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Extended Fourier analysis of signals 

Abstract. This summary of the doctoral thesis [8] is created to emphasize the close connection of the proposed 
spectral analysis method and the Discrete Fourier Transform (DFT), the most extensively studied and frequently 
used approach in the history of signal processing. It is shown that in typical application case, where the uniform 
data readings are transformed to the same number of uniformly spaced frequencies, the results of the classical DFT 
and proposed approach coincide. Performance differences appear if the length of the DFT is selected greater than 
the length of the data. DFT resolves the problem of unknown data by padding readings with zeros up to DFT length, 
while the proposed Extended DFT (EDFT) deals with this situation in a different way, it uses Fourier integral 
transform as a target and optimizes the basis for transformation in the extended frequency set without imposing 
restrictions on the time domain. Consequently, Inverse DFT (1DFT), which is suitable for EDFT results, gives not 
only known readings but also the extrapolated data where classical DFT can only be returned zeros, and higher 
resolution is reached at frequencies where the data has been successfully extrapolated. EDFT has been shown to 
able to process data with missing readings or gaps inside or even nonuniformly sampled data. Therefore, EDFT 
significantly extends the usability of DFT based methods, where previously these approaches have been considered 
as not applicable [ 10 - 45 ]. EDFT finds a solution in an iterative way that requires repeated calculations to obtain an 
adaptive basis, and that makes numerical complexity much higher compared with DFT. This disadvantage was a 
serious problem in the 1990s, when the method was proposed. Fortunately, since then, computer power has 
increased so much that the use of EDFT can be a real alternative nowadays. 


1 Introduction 


A Fourier transform is a powerful tool for signal analysis and representation of a real or complex¬ 
valued function of time x(t) (hereinafter referred to as the signal) in the frequency domain 

oo 

F(u>) = I x(t)e~j° >t dt, (1.1) 

j — OO 

x(t) = — [ F(co)e Jcot dco. (1.2) 

J — oo 

The Fourier transforms orthogonality property provide a basis for the signal selective frequency 
analysis 



e JW ° t el a)t dt — 2n8(co — u> 0 ), 


( 2 ) 


where a>, coo are cyclic frequencies, j is an imaginary number such that j 2 =- 1 and co-coo) is the 
Dirac delta function. Unfortunately, the Fourier transforms calculation according to (1.1) requiring 
knowledge of the signal x(t) as well as perfonning of integration operation in the infinite time 
interval. Therefore, for practical evaluation of (1.1) numerically, the observation period and the 
interval of integration is always limited by some finite value 0 and the signal is known in the time 
interval -0/2<t<0/2. The same applies to the Fourier analysis of the signal sampled versions - 
nonuniformly sampled signal x(4) or unifonnly sampled signal x(kT) for k=-co,.. .,-1,0,1,.. .,+qo. 
Only a finite length sequence x(fc) or x(kT), k= 0,1,2,. ,.,K- 1, are subject of Fourier analysis, where 
A is a discrete sequence length, T is sampling period, and the signal observation period is equal to 
0 =tK-i~to or 0=AT To avoid aliasing and satisfy the Nyquist limit, unifonn sampling of continuous 
time signals should be perfonned with the sampling period T<7dCl, where Q is the upper cyclic 
frequency of a signal x(t). Although nonuniform sampling has no such a strict limitation on the 
mean sampling period T S =Q/K, in the subsequent analysis we suppose that both sequences, x(&) 
and x{kT), are derived from a band-limited in Q signal x(t). Let's write the basic expressions of 
classical and extended Fourier analysis of continuous time signal x(t) and its sampled versions x(&) 
and x{kT). 
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2 Problem formulation 

“The formulation of a problem is often more essential than its solution which may be merely a 
matter of mathematical or experimental skill. To raise new questions, new’ possibilities, to regard 
old problems from a new angle requires creative imagination and marks real advances in science. ” 

Albert Einstein and Leopold Inf eld, Evolution of Physics, 1938. 

2.1 Basic expressions of classical Fourier analysis 

The classical Fourier analysis dealing with the following finite time Fourier transforms 


f 0/2 

F 0 (n>) = 1 x{t)e~^ t dt, 

0/2 

K-1 

(3.1a) 

F 0 (m) = V x(t k )e~ ja>tk , 

(3.1b) 

k= 0 

K-l 


F 0 (o>) = ^ x(kT)e~^ kT , 

(3.1c) 

k= 0 

1 f- Q 

*©(t) = 7 T- F @ (a))e ja>t d(jo, 

2n 

(3.2) 


where (3.2) is the inverse Fourier transform obtained from (1.2) for a band-limited in Q signal. 
Transforms (3.1b) and (3.1c) are kn own as Discrete Time Fourier Transforms (DTFT) of the 
nonuniformly and uniformly sampled signals. The reconstructed signal xe(t) outside the 
observation period 0 vanishes quickly reaching values close to zeros. The signal amplitude 
spectrum is the Fourier transfonn (3.1) divided by the observation period 

S 0 (m)=^F 0 (m). (4) 

The frequency resolution of the classical Fourier analysis is inversely proportional to the 
observation period 0, thus, the longer interval of signal analysis, the higher resolution is achieved. 
Obviously, one can get the fonnula (3.1a) by truncation of infinite integration limits in (1.1) and 
the DTFT (3.1a) and (3.1b) in a result of replacement of infinite sums by finite ones. This mean, 
the classical Fourier analysis supposed that the signal outside 0 is zeros. In other words, the Fourier 
transfonn calculation by fonnulas (3.1) is well justified if applied to time-limited within 0 signals. 
On the other hand, a band-limited in Q signal cannot be also time-limited and obviously have 
nonzero values outside 0. Generally, the Fourier analysis results obtained by using the exponential 
basis tend to the Fourier transform, if 0— »oo, while in any finite 0 there may exist another transform 
basis providing a more accurate estimation of (1.1). 

2.2 Basic expressions of extended Fourier analysis 

The idea of extended Fourier analysis is finding the transfonn basis, applicable to a band-limited 
signals registered in the finite time interval 0 and providing the results as close as possible in terms 
of the Z. 2 -norm (or the Euclidean norm) to the Fourier transfonn (1.1) defined in the infinite time 
interval. The fonnulas for proposed extended Fourier analysis could be written as 



r 0/2 


F a (<u) = 

x(t)a(n>, t)dt, 

'-0/2 

K-1 

(5.1a) 

F a (<u) = 

^ x(t fc )a(m, t k ), 

k= 0 

(5.1b) 
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K-l 

F a (oo) = ^ x(kT)a(co,kT ), (5.1c) 

k =o 

1 f n 

x a (t) = — \ F a (( 0 )e J0)t d( 0 , (5.2) 

2n J_ n 

where in general case the transform basis a( co,t), a(co,tk) and oi co,kT) are not equal to the classical 
ones (3.1). Note that the inverse Fourier transfonn (5.2) still preserves the exponential basis and 

the Parseval-Plancherel formula j_fx a (t)\ 2 dt — ^ J^|F a (o))| 2 dco holds for it. 

To ensure that the results of transforms (5.1) are close to the result of the Fourier transform (1.1) 
for the signal x{t), the following minimum least squares expression will be composed and solved 

|F(u>) — F a (u>) | 2 -* min. (6) 

Unfortunately, as already stated above, the calculation of F(cd) cannot be performed directly for a 
band-limited signal. So, to compose (6), we should find an adequate substitution. Let's recall that a 
complex exponent at cyclic frequency coo and with a complex amplitude S( ox) is defined in the 
infinite time interval as 


x(co 0 ,t ) = S(u> 0 )e- ,w ° t , —oo < t < oo. 

The Fourier transfonn of a signal (7) can be expressed by the Dirac delta function (2) 


L 


x(u> 0 , t)e J(x>t dt = 2nS(co 0 )S(co — o> 0 ). 


(V) 

( 8 ) 


Now, we will use (7) as a signal model with known amplitude spectrum S(coo) for frequencies in 
the range -Q.<coo<Q and, in the expression (6), substitute F{cd) by the Fourier transfonn of the 
signal model (8) and the signals x{t), x(4) and x(kT) in (5.1) by the signal models (7), 
correspondingly. Finally, the integral least square error estimators for all the three signal cases get 
the fonn 


n 


n 






27T S(cl> 0 )S(co 


Q 


Q 




-n 


27tS(co 0 )S(co 



r®/ 2 

2 


M 0 ) - 

1 S((jo 0 )e J(I>ot a((jL),t)dt 

cI(a)q, 

(9a) 


7-0/2 




K-l 

2 


- co 0 ) - 

- V S((n 0 )e ]a>otk a((n, t k ) 

d(j)Q, 

(9b) 


k =0 




K-l 

2 


Mo) - 

5(o) 0 )e- ,w ° ,cT ci:(a), kT) 

d (0 q . 

(9c) 


k =0 


The solutions of (9) for a definite signal model (7) provide the basis a( o,t), a{o),tk ) and a( ooJiT) for 
the extended Fourier transforms (5.1). To control how close the selected signal model amplitudes 
S{ ox) are to the signals x{t), x(4) and x(kT) amplitude spectrum, we will find the fonnulas for 
estimate signal amplitude spectrum S f/ (oo) in the extended Fourier basis a(co,t), a{co,tk ) and a(<x>,kT). 
The fonnula (8) is showing the connection between the signal model Fourier transfonn and its 
amplitude spectrum, from where Sica)) could be expressed as signal model Fourier transfonn 
divided by Ind^co-coo). Taking (8) into account, S„{o)) is calculated as the transforms (5.1) divided 
by the estimate of InS^oo-coo) in the extended Fourier basis, which is detemiined from (9) in the 
case A=0 and ox)=o), 

r ©/2 


F a (^) = 


F a (m) = 


f-e/2 eiaJt( x{oj, t)dt 

Xfc = QX(t fc )g(6J, t fc ) 

Zfc=o e jatk a{(0, t k y 


(10a) 


(10b) 
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_ Zk=o x(kT)a((o,kT) 
eja>kT a((o, kT)’ 


(10c) 


and showing that the amplitude spectrum on the frequency co is estimated as ratio of the signal 
extended Fourier transform to the transfonn of exponent with a unit amplitude in the same basis. 
This is true also for the classical Fourier transform. For example, after substituting exponential 
basis a(co, t ) = e~ J0>t in (10a), the denominator becomes equal to 0 as in formula (4) for the 
classical Fourier analysis. 

The values of the denominator in formulas (10) are in inverse ratio to the frequency resolution of 
the extended Fourier transfonn. 

Before finding the extended basis functions for arbitrary S(coo), it is reasonable to consider a simple 
signal model having a rectangular form, S(axi) = I for -Q<oo<Q and zeros outside. Then the 
estimators (9) reduce to 


t 

-n 


rQ/2 

2 


^ j 


2n8((ji) — <u 0 ) — 

1 e JC ° ot a((L),t)dt 

d(i) Q; 

(1 la) 

J 

-Q 


2-0/2 




r n 


K- 1 

r—i 

2 


A= 

f 

2n8((jL) - ui 0 ) - 

> e^° tk a((jL>,t k ) 

d(j) Q; 

(lib) 


'-n 


L—l 

k= 0 




f Q 


K-l 

v—1 

2 


a =j 

1 

-n 

2n8((ji) — o> 0 ) — 

y e 2«o kT a(oo, kT) 

is — n 

d (0 q . 

(11c) 


The solution (11) allows us to establish a relationship between the classical and extended Fourier 
transforms. 


3 Problem solution 

In this section the integral least square error estimators (9) and (11) are solved and subsequent 
analysis of the obtained results is perfonned to find out only those solutions that can lead to 
practically realizable algorithms. 


3.1 Extended Fourier transform of continuous time signals 

The solution of (11a) for continuous time signal x(t) is found as a partial derivation 
= 0 , — j < t < j, and leads to the linear integral equation 


9 a(o>,T) 


[ 


2 ' 
0/2 


sin(Q(t — t)) 


a(cL>,t)dt — e i 0JT . 


( 12 ) 


'-0/2 71 O' ~ T ) 

Step by step solution of (12) is given in [4], Finally, the basis a(co,t) are obtained by applying a 
specific function system - a prolate spheroidal wave functions [1] y/k(t), k= 0,1,2,..., and are written 
as series expansion 

oo 

B k (<w) ... 

^(0- 


a(co, t) = ^ ; 

k=0 


Ah 


(13) 


The extended Fourier Transfonn of continuous time signal x(t) are given by 

OO 

E a (<u) = ^ B k (o))a k , -Q < a < 0, 


k =0 


Xa(t) = ^ W k (t)a k ’ -oo < t < OO, 


(14.1) 

(14.2) 


k =0 
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s a (<w) = 


Y'k=oB k (w) a k 

S” =0 |5/c(^)l 2 ' 


(14.3) 


where a k = X k = if%(t)dt, B fc (<u) = ^ k (m^)(-;) /c and 


- 0/2 


,L fe 0 


2fT 


the Parseval-Plancherel equality gives /^Jx a (t) | 2 dt = ^/^|F a (<n)| 2 dm = £“ =0 |a fe | 2 . 


The extended Fourier transfonn in accordance with (14.1) requesting a calculation of infinite sums, 
this mean, an infinite quantity of mathematical operations, therefore it’s impossible for real world 
applications. Theoretically, the value of the denominator £fc=olF/c (n>) I 2 in the amplitude spectrum 
formula (14.3) tends to infinite for K—> oo, and the extended Fourier transfonn (14.1) provide a 
super-resolution - an ability to determine the Fourier transfonn for the sum of sinusoids or complex 
exponents, if frequencies of them differ by arbitrary small finite value. 


3.2 Extended Discrete Time Fourier Transform 

In this subsection the minimum least square enor estimators (9b,c) and (llb,c) are solved and the 
extended Fourier transforms for uniformly and nonunifonnly sampled complex-valued signals are 
obtained. The proposed approaches have been developed in articles [5] and [6], where the 
derivations for real-valued discrete signals are given. 

The following notations are used in the matrix equations: superscripts X" 1 , X 7 , X* and X H denote 
inverse, transpose, complex conjugate and complex conjugate (Hennitian) transpose of the matrix 
X; ./ represents element-by-element division of two matrices with the same size; sum(X ) means 
addition of all matrix X elements and the diag(X ) forms the row vector by extracting the main 
diagonal elements from quadratic matrix X or it puts the elements of vector X on the main diagonal 
to form a diagonal matrix. 


3.2.1 Particular solution for discrete time signals 

The solutions of (llb,c) can be obtained similarly to (11a), as partial derivatives of 


5A 


and 


aA 


da(u>,ti) 


= 0 


da(cj,lT ) 


= 0 for /=0,1,2,...,/f-1, and leads to the systems of linear equations 


K-l 

I 

k =0 
K-l 

I 

k =0 


sin(fl(t fc - t t )) 
n(t k - tO 

sin(Q(/c — l)T) 

n(k - Tyf 


a(co, t k ) = e 


a((jo,kT ) = e - ,wiT . 


The solution of (15) in the matrix fonn is expressed as 

A = R -1 F 


(15a) 


(15b) 


(16) 


where \o>(Kx 1) and E®(Xxl) are the extended Fourier and the exponential basis. 

The formulas of Extended Discrete Time Fourier Transfonn (EDTFT) for signal model <S(oo) = l, 
-Q<cm<Q are derived by substituting of transformation basis (16) into expressions (5) and (10) 

F a (co) — xR -1 E w , —Q < 6J < _Q, (17.1) 

x a(0 — xR -1 E t , —oo < t < oo, (17.2) 

xR _1 E £l) 

S a (oj) = — —(17.3) 

aV J pUD-lp V 7 

The matrices for nonunifonnly sampled signal x(tt) are composed as follows 
x(lxX): x(&), E« (Kxl): e~^\ R (KxK): r ljk = Sin( ^ fe ~ t;)) , E t (Kx 1): 

Uniformly sampled sequence x{kT) could be considered as a special case of nonuniform 
sampling at time moments tic=kT, k=0, 11, then the matrices in (16, 17) are fonned as 
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x(l xK): x(kT), E® (Kxl): e~ joilT , R ( KxK ): r lk = 


sin(Q(k-0T) 


, E,(Axl): e, = 


sin (Q(t-IT)) 


n(k-l)T ’ v ' 1 n(t—lT) 

If sampling of signal x(kT) is done with Nyquist rate, T-niil, the matrix R becomes a unit matrix I 
and the fonnula (17.1) coincide with classical DTFT (3.1c), but the fonnula (17.3) reduces to the 
well-known relationship between discrete signal Fourier transform and its amplitude spectrum 

F a (m)=F 0 (m)=x E w> (18.1) 

= ^XE W . ( 18 - 2 ) 


Whereas for nonunifonnly sampled signal x(&-) the matrix R^I, even if mean sampling period 
T s =nlD- and formulas (17) give the results that are close to unifonn case and superior to those 
obtained by the classical nonunifomi DTFT (3.1b). The resolution by frequency in both sampling 
cases equals to 1/FT, which is a nonnal frequency resolution. While for oversampled signals, T (or 
T s ) < /z/Q, the EDTFT approach can provide a high frequency resolution and improved spectral 
estimation quality. Unfortunate an achievement of such results is limited by finite precision in the 
mathematical calculations and by restrictions on frequency range in the process of signal sampling. 
The theoretical value of the denominator in (17.3) E^R _1 E W = K and the frequency resolution 
should increase proportionally to the number of samples in the signal observation period 0. In the 
border-case, if the number of samples within 0 is increasing infinity, A—>• oo, and the discrete time 
signal tends to the continuous time signal x(t), the EDTFT (17.1) gives the same result as (14.1). 


3.2.2 Generalized solution for discrete time signals 

Now, we will consider the solution of the minimum least square error estimators (9b,c) for 
arbitrary selected signal model S(coo). The derivation fonnulas for both estimators are like the 
ones given in the previous section. For example, a partial derivation of (9b) by the basis, 


9A 


9 


= 0 for 1=0, 11, provides the least square solution 


j ^27tS(<x> 0 )<5(u> — u> 0 ) — ^ 5(o) 0 )e- ,w ° tfe ct:(a), 5*(o) 0 )e - ,w ° ti d(x> 0 = 0, (19) 


Equation (19) can be rewritten as 
' -n 


^ fj lS(co 0 )l 2 e JOJ ° ( ' tk ~ t ^da> 0 ^j a(a>, t k ) — 27T J \S(co 0 )\ 2 e~^ a>otl 8((jL> — (jL> 0 )dco 0 . (20) 

k= o' ~ n ' 

The filtering feature of the Dirac delta function f f(x)8(x — x 0 )dx — f(x 0 ) applied to the 

right part of (20) gives the final form of the system of linear equations for /=0,l,2,...,FiT, 
k -i n 

Xl^j \ s M\ 2 e J “° (tk - tl) doo 0 ja(oo,t k ) = \S(co)\ 2 e~j^, (21a) 

k=0 

K - 1 / Q . 

Xl^j l s M\ 2 e JcO0(k -° T dco 0 )a(a>,t k ) = \S((o)\ 2 e~^ lT , (21b) 

k = 0 

where |5(<u)| 2 is the signal model power at ax=ax The system of linear equations (21b) is 
applicable for unifonnly sampled signal x(kT) and can be derived from (9c) in a similar way as 
(21a). 

The EDTFT basis a(a),tk) or o{oxkT) can be found by applying different solution algorithms to the 
system of linear equations (21). In general, basis A® (Kx 1) is obtained in the matrix fonn as 

A w = iSCnOfR-X (22) 

and inserting (22) into expressions (5) and (10) yields the fonnulas for calculation of the EDTFT 
F a (co) — xA w = |5 (<u)| 2 xR _ 1 E w , —Q < ox < Q, (23.1) 
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x a (t) — X R — 00 < t < oo, 


S a (oj) = 


xA w x|5(o))| 2 R 1 E 6 


F H A 


ESIS(m)|2R-iE 6 


x R~ 1 e 6J 

EgR"^*,' 


(23.2) 

(23.3) 


The elements of the matrices R (KxK) and E, (Kx 1) in the formulas (22, 23) are expressed by 
integrals 

n,k = J~f |S(w 0 )|V w ° (tk_t!) dm 0 or r Uk = [ |S(m 0 )|V w ° (k -° T dm 0 , (24.1) 

2n J_ n 2n J_ a 


e t —— [ "|5(m)| 2 e- ,w ^ tl) dco or e l —— [ \S((jj)\ 2 e^^ lT ^doj, (24.2) 

2nJ_ Q 2n J_ n 

for nonunifonnly or unifonnly sampled signal cases, respectively. If the signal and its model power 
spectra are close, |5’ a (<n 0 )| 2 * |S(u> 0 )| 2 , then (24.1) is also an estimate of the autocorrelation 
function of the sequence x. The inverse transfonn (23.2) calculated on time moments t=tkO r t=kT, 
k=0,l,2,...,K-l, returns back the input sequence x undistorted, as the elements of matrices E, 
become equal to R. Case signal model S(coo)=l the fonnulas (22) and (23) reduce to (16) and (17). 
The frequency resolution of the EDTFT is in inverse ration to |S(u>)| 2 E^R -1 E w and varied in 
the frequency range -Q<<u<Q. 


3.3.3 Iterative EDTFT algorithm 

Calculation of the EDTFT by fonnulas (23) requires knowledge of the signal model spectrum 
which generally is not known. At the same time, the amplitude spectrum obtained in the previous 
section by the fonnula (17.3) can be used as a source of such information. This suggests the 
following iterative algorithm introduced in [5], where the signal model spectrum S(coo) tends to the 
signal spectrum S a (co): 

Iteration 1: Calculate S a (oj) (17.3) applying default signal model - 5 '(od)=1. 

Iteration 2: Calculate S^ 2 \co) (23.3) by using the signal model - 5^(<u 0 ). 

f3^ (2) 

Iteration 3: Calculate S a (oj) (23.3) by using the signal model - S a (u> 0 ). 


Iteration I: Calculate (oj) (23.3) by using the signal model - ±J (<A))- 

The iterations are repeated until the given maximum iteration number is reached or the power 


(f-i), 


spectrum does not alter from iteration to iteration - (oj) | 


Ti-i) 


(a.) | . 


The EDTFT output F a (co) (23.1) is calculated for the last performed iteration I. 

By default, the signal model S( «>)= 1 is used as input for the EDTFT algorithm. However, 
additional in formation about the signal to be analyzed can be applied to create a more realistic 
signal model for the EDTFT input and to reduce the number of iterations required to reach the 
stopping iteration criteria. 


4 Extended DFT 

EDTFT considered in the previous section is a function of the continuous frequency (-Q<a<Q), 
while describing below EDFT algorithm calculate EDTFT on a discrete frequency 
set -Q< 00 n<O. for n= 0,1,2,... ,N- 1. The number of frequency points N>K and it should be selected 
sufficiently great to substitute the integrals (24.1) used for calculation of the matrix R ( KxK) in the 
expressions (22, 23) by the finite sums 
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_Q n=0 

fi N-i 

fi,/c = ^:J l^(^o)l 2 e ; " o(k “° T rf^o ~ 


2 pj^nO'k tl) 


2gjton(k QT 


IJf=0 9 \,2 9 ... 9 K-\. The matrices composed of (25.1) and (25.2), 

r o,o(0) r o,i( f i — ^o) r o,*:-i(t*:-i — to) 

r i,o(to — ti) rvi(O) t i,k-i(Pk-i ~ t]_) 


R = 


R = 


- r *:-i,o(to — tn-i) r K-i,i(.ti — t K _ 1 ) 
r o,o(0) r 0f iCO 

r i,o( — T) r lf i(0) 


r x-i,x-i(0) 
r o,K-i(.(.K ~ 1)0 
r^-iCC^ - 2)T) 


(25.1) 

(25.2) 


(26.1) 


(26.2) 


L^at— i,o( — 1)T) r K-i,i ( — (K — 2 )T) ■■■ r K-i,K-i(0) 

possess Hermitian symmetry, r l k — r@, but (26.2) for a unifonnly sampled signal has also a 
Toeplitz structure. The matrix elements r i k represents the autocorrelation function and can be 
calculated by applying the IDFT to the signal model power spectrum |S(co n )| 2 . The frequency 
Q/^=2/ ) ',=/\ in (25), where/, is the signal upper frequency and /v is the Nyquist rate of a band- 
limited signal, and it is assumed to be normalized (equal to 1) in DFT calculations. The choice of 
the frequencies {<x>n}={ 2;/,} depends on the number of frequencies needed for accurate estimation 
of (25) as well as for detailed signal spectrum representation, and the limitations on the total amount 
of computation. Eventually, the uniform set of frequencies in range [-Q;Q] is preferable in most 
application cases. 

The EDFT may be expressed by the iterative algorithm 


R6) = 1eW®E h , 

N 

(27.1) 

p(0 = xA (0 = x(r6)) -1 EW6), 

(27.2) 

, m _ *(r“T 1 e. 

(27.3) 

diag(E H 


W6 +1 ) = diag (\S^\ 2 ^j, 

(27.4) 


for iteration number z'=l ,2,3,...,/, wherein (27.1) is the sum (25) in matrix fonn. The matrix E (KxN) 
has elements e~ ]2n f ntk or e~ j2n ^ nkT casQ sampling of x done uniformly. By default, the diagonal 
weight matrix W (,) (NxN) for the first iteration is a unit matrix, W (I) =I. If the other diagonal matrix 
is used as input for the EDFT algorithm, it should have at least K non-zero elements. For the next 
iterations W (,+1 f is filled with power spectrum values calculated by (27.4). There could be additional 
criteria for stopping the iterations before the maximum number of iterations / is reached, for 
example, the iterations could be interrupted, if the relative change in the power spectrum, 
\sum(W {i+l) )-sum(W (i) )\/sum(W {2) ) for i> 1, is smaller than a given threshold. 

The IDFT may be applied to output F of each iteration and returns original K samples of unifonn 
or nonunifonn sequence 

x = 2 fe / (28) 


Since the length of the frequency set N>K, then (28) can be modified to obtain an extrapolated 
sequence x a (1 x/V) - x a (t m ) ovxJjnT), m=0,l,2,...JV-\, 


x 


a ~ 



H 

N> 


(29) 
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where exponents matrix E« (NxN) has elements e~ j2n ^ ntm or e~ ]2n ^ nmT case of uniform x a , and 
xx H < x a Xa — FF /y according to Parseval-Plancherel theorem. Reconstructed by the formula 
(29) sequence is the original sequence plus forward and backward extrapolation of x to length N 
and/or interpolation if there are gaps inside of x. The maximum frequency resolution is limited by 
the length N of frequency set, not by the length K of sequence x as in the application of classical 
DFT. It means, the EDFT can increase the frequency resolution NIK times in comparison with the 
classical DFT. This can be verified by comparing the diagonal elements of the product of IDFT and 

DFT basis, diag (j^E h E^ , which are equal to KIN at all frequencies, with the relationship, 

0 < diag Q-E h A^ = ~F./S < 1, corresponding to the IDFT and EDFT basis A (27.2). At the 

same time there is a restriction on the frequency resolution sum(¥ ,/S)=NK, which is satisfied by 
iteration, and to achieve a high resolution at certain frequencies, the EDFT must decrease the 
resolution on other frequencies. The deviation \sum(F./S)-NK\ also could be used as an additional 
criterion for stopping of iterations, because indicates the possible inaccuracy in the obtained results, 
mainly caused by the finite precision in calculations. If this happens, the result of the previous 
EDFT iteration should be considered as a final one. 

In a border-case N=K, the iterative algorithm output does not depend on weight matrix W and 
the optimal EDFT basis is found in a non-iterative way (in a result of the first iteration) [7]. 

5 EDFT and other nonparametric approaches 

In the previous sections, starting with the Fourier integral (1) and using its orthogonality property 
(2), by establishing and solving the minimum least square error estimators (9), the Extended DFT 
is obtained analytically. In the next, a comparison with known nonparametric approaches - 
Capon filter, Generalized (Weighted) Least Squares (GWLS) solution and High-Resolution 
Discrete Fourier Transform (HRDFT) introduced by Sacchi, Ulrych and Walker in 1998 are 
done, and the ways and opportunities of derivation of an iterative EDFT algorithm based on 
these methods are analyzed briefly. 

5.1 Capon filter approach 

The Capon filter kn own also as Minimum Variance spectral estimate (see [3,10,11, 24]) can be 
viewed as the output of a ha nk of filters with each filter centered at one of the analysis 
frequencies 

y w OfO = Y£= oM( n ~ k')T)h (a QcT') = xh w , n = 0,1,2,... . (30) 

In the matrix notationx = [x(nT), x((n — 1)T), ...,x(fn — K + l)T)j is the filter input signal 

and h w = [/i w (0),h w (70. h^K - 1 )T)] t is the filter coefficients. Here the subscript co 

indicate a dependence on the filter’s center frequency. 

The Capon filter is designed to minimize the variance on the filter output 
= £{ly w OfOI 2 } = £{y"(nr)y 6) (n7’)} = E{h^x H xhJ 

= hMx H x}h w = h"R,h w , P ; 

subject to the constraint that its frequency response at the frequency of interest co has unity gain 

K-l 


H(co) = 

V h^kT)e-^ kT = E£h w 

= 1 , 

(32.1) 


k =0 

K-l 



H(a)) = 

Y hl(kT)ei“" = h!X 

= 1 , 

(32.2) 


k =0 




Where e{.} denotes the expectation operator and the matrix E® (Kx 1) has elements e . The 
constraints (32.1) and (32.2) must be satisfied by filter (30) and Hennitian transpose filter 
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y% ( nT ) = h^x H , correspondingly. The matrix R = f{x H x} (KxK) is the sample autocorrelation 
matrix and it can be composed of the values of the signal autocorrelation function. For example, so 
called biased estimate is calculated by 


K-l-l 


r xx m = \ ^ x((k + l)T) x*{kT ), Z = 0,1,2. K- 1, 

k=0 

and, considering that r xx {—lT) — r xx (lT ), the sample autocorrelation matrix is filled as 
r o,o(0) r-oA-T) ___ r 0 «•_,(—(if — 1)T)] 

R _ r i,o(T) r i,i(0) ^ 1 ,/c—i( (Zf — 2)T) 


(33) 


(34) 


k-i,o((/f - l)n r K _ 1 x ((/f — 2)T) ••• t"x-i,x-i(0) 

Mathematically, the Capon filter coefficients can be obtained by minimizing the variance (31) under 
the constraints given by (32.1) and (32.2) 

J = h£R*h w - g(ElK - 1 ) - MKK - 1 ) -> min, (35) 


where p,X are Lagrange multipliers. The conditions — 0 and —77 = 0 must be fulfilled to 

on,. 


dK, 


detennine the minimum of (35). Both requirements lead to the same solution 


R@E 


h<D 


F T R _1 F* 


and, traditionally, the Capon power spectrum is computed as 

u - 1 

Pcaponi.^') ~ h, 


t D U — 

l co 


E T R _1 E*' 


(36) 


(37) 


To obtain an iterative EDFT algorithm from the original Capon filter approach, the sample 
autocorrelation matrix R v (34) must be substituted by R /= E WE 7 . The matrix R 7 ( KxK) can also 
be obtained as a transpose of the EDFT matrix R defined by (26). The elements of quadratic 
diagonal matrix W (NxN) represent an estimate of power at time moment nT= 0, detennined from 
one sample at output of each Capon filter 

x(R t ) 1 E^ 


|y w ( 0)| 2 = |xhj 2 = 


eur t )~ 1 e: 


(38) 


- l OJ J CO 

where the filter input sequence x (30) is related to the EDFT input sequence x 
x(nT) — x((K + k — 1)T) or x(t k ) — x(t K+k _ 1 ), k=0,-l,-2,..,-(K-l), for unifonnly 
nonunifonnly sampled sequence cases, respectively. 

Eventually, an iterative algorithm can be fonned as follows 


as 

or 


R r © = E*W®E r , 

(39.1) 

„(0 x(R r <')) _1 E*. 

(39.2) 

Cavon diag( E 7 ’(R 7 ’©)-i E *)' 

w(i+1) = diag ( Sf ap0n Y 

(39.3) 


with the initial condition for W (1) =I and the iteration number i=l, 2,3,.../. The estimate of the power 
2 


spectrum 


5 (0 

*Capon 


different. It shou 


coincides with the results of the EDFT, while the phase spectrum, definitely, is 

d be noted that the calculation of the Capon filter output power by (37) is 
theoretically well justified, whereas the derivation of (39) requires ad hoc assumptions and 
substitutions, and actually is a measurement of power obtained from just a one sample at the output 
of the filter. This leads to conclusion that the approach (39) is simply a filter-bank interpretation of 
the EDFT, similarly to the DFT which can also be considered as a bank of filters. In addition, an 
iterative algorithm derived based on Capon filter cannot reveal all the EDFT capacity, such as the 
ability to estimate DFT (27.2) and restore the signal (28, 29). 
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5.2 GWLS solution 


The Generalized (Weighted) Least Squares approach (see [3,15,18, 34]) in the spectral analysis 
could be based on the following data model 

x r — E^S GVVii s((x>) + e^, (40) 

with e <2 denoting the noise and interference (signals at frequencies other than a>) component, 
and E l)S GWLS (o)) representing the signal component on the frequency of interest with unknown 
complex amplitude S gwls (oj). The GWLS minimizes 

[x T — E^S’ GM / LS ((d)] h Q W -E'vSgwlsM], (41) 


which is solved by 


J GWLS 


(«) = 


e^Q 


-1 


e LQ _1 ej 


(42) 


where Q (KxK) is the covariance matrix of the data model component e Q . There are two special 
cases of GWLS called Weighted Least Squares (WLS) and ordinary Least Squares (LS). WLS 
occur when all the off-diagonal entries of Q are 0, while LS solution is obtained from the GWLS 
under the assumption that e Q at (40) is a white noise, hence Q=I. 

The problem of GWLS estimator is that, in general, the noise covariance matrix Q is not known, 
and must be estimated from the data along with the Sgwls(co). The initial estimate (the 1 st 
iteration) could be equal to LS solution, it is (42) with Q=I. Next, to ensure that the GWLS 
solution works in an iterative way as EDFT do, the covariance matrix should be calculated as 
Q=R 7 =E*WE r under the assumption W = diag(\S GW ls (oj)\ 2 ). In result, GWLS solution (42) 
coincides with the EDTFT formula (23.3) 


El ) (R T )- 1 x r xR -1 E 6 


El ) (R 7 ’)~ 1 E* 


pUD-lp 


= S a ((o). 


(43) 


$gwls (aO — 

L (x) J ‘-‘CO 

and, as shown in the Section 3.3.3, can be successfully used to update of the amplitude spectrum 
iteratively. 

Although substitution of a noise covariance matrix by R r would be easy done, it is not supported 
by GWLS data model (40), from where the matrix Q represents the data model component e Q 
only and the signal component E^S GlViS (n>) must be excluded from it, whereas the matrix R 7 
is calculated for the entire signal x T , including e Q and E^S GWLS (a>) . Furthermore, the 
derivation of EDFT in the Section 4 shows that the signal is restored by applying IDFT (28) to 
the Extended Fourier transform, x T — ^ E * F A E*S, and not as an inverse of the Amplitude 


spectrum (27.3), which is a scaled version of (27.2) by a frequency dependent weight factor. Using 
an estimate SGWLs(cd) = Sa(cd) in the data model (40) leads to a predetermined split of the signal at 
frequency co in between both components, where the noise part may be expressed as a difference 
of EDFT outputs F and S. The conclusion reached is that making the derivation of the Extended 
DFT algorithm possible, invalidates GWLS minimization expression (41) which require 
separation of both data model components. 


5.3 High-Resolution DFT 

The third method considered here is High-Resolution DFT proposed by Sacchi, Ulrych and 
Walker in [9]. The authors presented an iterative nonparametric approach of spectral estimation, 
which minimizes the cost function deduced from Bayes’ theorem and, as well as Extended DFT, 
makes it possible to obtain high-resolution Fourier spectrum. The HRDFT algorithm can be 
reduced to the following iterative procedure: 

r(0 = Iew^E h , (44.1) 

F f RDFT = x(R®) _1 EW© (44.2) 
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W/( l+1 ') — diag 


_ F (0 

r HRDFT 


(44.3) 


for iteration number z'=l,2,3,.../and with the initial condition YV (I, =I. 

The IDFT (28) applied to any iteration output (44.2) returns back the sequence x undistorted. The 
main difference between approaches is that the HRDFT algorithm lack of formula to estimate of 
amplitude spectrum (27.3). Instead, as input for the next iteration, it uses the Fourier spectrum 
estimated in the previous iteration. Thus, the results of FIRDFT differ from output of EDFT 
significantly. FIRDFT iterates to the solution where the signal is approximated by K frequencies 
while the power on other N-K frequencies becomes negligible. Each valuable frequency is resolved 
with maximum resolution restricted by the length of HRDFT. Also, it still obeys the same limit on 
the sum of resolutions by frequency ( KN) as DFT and EDFT. 

The authors [39] investigated algorithms having different from (44.3) weights for adaptation of the 
correlation matrix (44.1), although only the amplitude spectrum (10) derived accordingly to the 
minimum least squares expression (9) and calculated by (27.4) fits perfectly to an iterative update 
of the matrix R and returns results that are closest to the Fourier transfonn in the Z 2 -nonn sense. 


6 Computer simulations 

The EDFT algorithm is validated on the data which are similar to those that have been used in [5, 
7, 8]. The true spectrum of the first test signal consists of a band-limited noise (flat) in the frequency 
range [-0.5...-0.25] Hz, a rectangular pulse in the range [0...0.25] Hz and two unit power complex 
exponents at frequencies 0.35 Hz and 0.3985 Hz. These three components represent random, 
transient pulse and deterministic parts of a composite signal with the upper frequency@=0.5 Hz. 
Unifonn and nonunifonn sequences of the length K= 64 samples are derived by simulating 10-bit 
Analog-to-Digital Converter (ADC). Sampling and mean sampling periods of both sequences are 
equal to 1 second, T=T S = Is. Sampling time points for the nonunifonn sequence are generated as, 
tk=kT+zk, k= 0,1,2,...@-1, where {»} are unifonnly distributed random values in the range 
[0...0.8s]. Thus, the true spectrum of sequences consists of three non-overlapping in frequency 
domain components and ADC added floor noise (~-60dB), and it is symbolized by red color lines 
in the Figures 1-4. 

The plots in Figures 1 and 2 show the perfonnance of EDFT (black lines) for uniform and 
nonuniform sequences and allows to compare it with the classical DFT (blue lines). The number of 
frequencies (the length of DFT) is chosen equal to N= 1000, which gives spectral estimates with 
DFT frequency bin spacing 2fJN=0.00 1 Hz. This means that the range [-0.5...0.5[ Hz is unifonnly 
covered by frequencies and used for the calculations in (25, 27) and for the signal representation in 
the frequency domain (spectral plots). Figures la and 2a display the power spectra of EDFT 
calculated as 101og(|S| 2 ) in a non-iterative way. The input matrix W in this case is composed of 
values of the true spectrum (red line in the plots), therefore there is no need for further iterations. 
Non-iterative estimate is very close to the EDFT 15 th iteration depicted in Figures 1 b and 2b, where 
the matrix W=I used in the input and confirms the correctness of the iterative algorithm. The Figure 
lc (2c) shows the Power Spectral Density (PSD) calculated by the EDFT as 101og(|F| 2 /A) and 
proves the expectations, that PSD estimate on a complex exponent should increase in a value in 
comparison with the classical DFT if the proposed method achieves a high resolution around this 
frequency. 

Figure Id and 2d plotting the relative frequency resolution for the EDFT 15 th iteration calculated 

as —-— F./S (Id) or —-— F./S (2d) in respect to the DFT for which, in accordance with (18), it is 
2f u TK 2 At FA 

simply equal to 1 at all frequencies. The value 2f u T=2f u T s =\ and this means that the signal is 
processed in one Nyquist zone. The DFT is showing a normal frequency resolution, whereas the 
EDFT have ability to increase the resolution (in plot appears values >1) around the powerful signal 
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components and decrease the resolution (in plot appears values <1) at frequencies where the signal 
has weak power components. 

(a) 



(b) 




(d) 



Figure 1. Uniform complex-value sequence - the estimate of: 

(a) Power spectrum - True (red), DFT (blue) and non-iterative EDFT (black), 

(b) Power spectrum - True (red), DFT (blue) and EDFT (15 th iteration, black), 

(c) Power Spectral Density - True (red), DFT (blue) and EDFT (15 th iteration, black), 

(d) Relative frequency resolution - DFT (blue) and EDFT (15 th iteration, black). 


The EDFT is called as a high-resolution method and that is true, but with the following remark - it 
keeps the same 'summary' resolution as the traditional DFT or, in other words, squares under black 
and blue curves in Figure Id (2d) are equal. The maximum frequency resolution is limited by value 
of division NIK. For example, if K=6A and N= 1000, then the EDFT can potentially improve the 
frequency resolution 1000/64«16 times. The peak resolution is achieved on a deterministic signal 
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part - at frequency 0.35 Hz. The resolution on 0.3985 Hz exponent does not reach maximum value 
because of its frequency is not on EDFT grid (0.001 Hz) and the estimate splits between two 
adjacent frequency bins, which is known artifact of the DFT analysis. 

(a) 




(c) 




Figure 2. Nonunifonn complex-value sequence - the estimate of: 

(a) Power spectrum - True (red), DFT (blue) and non-iterative EDFT (black), 

(b) Power spectrum - True (red), DFT (blue) and EDFT (15 th iteration, black), 

(c) Power Spectral Density - True (red), DFT (blue) and EDFT (15 th iteration, black), 

(d) Relative frequency resolution - DFT (blue) and EDFT (15 th iteration, black). 

Pulse signal ([0...0.25] Hz) is processed by EDFT with about the same resolution as DFT (~1). The 
relative resolution of random component ([-0.5...-0.25] Hz) fluctuates around 1, while in the 
regions where just ADC noise can be found, EDFT decreases the frequency resolution bellow the 
normal. The simulation showed that EDFT can successfully estimate random, transient and 
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deterministic signal spectra and provide results superior to those produced by traditional DFT. 
EDFT estimates in Figures 1 and 2 are close to each other and proves that the proposed approach 
can handle unifonn and nonuniform sequences with the same quality, while the efficiency of 
classical DFT gets worse in case of nonunifonn data. 

Figure 3 explains the difference in performance between unifonn and nonunifonn inputs, where 
the spectra of unifonn and nonunifonn sequences are analyzed in the extended frequency range, 
[-l...l[ Hz. The number of frequency points and the upper frequency are increased two times, 
,¥=2000 and @=1 Hz. This means that the step by frequency remains the same as in the previous 
plots. The true spectrum of sequences at frequencies above 0.5 Hz consists only of floor noise 
(«-60dB) added by ADC. The actual result depicted in Figure 3a shows periodicity of the DFT and 
EDFT spectral estimates, which cannot be avoided for unifonn sequences. In contrast, EDFT 
applied to the nonunifonn sequence returns conect power spectrum in Figure 3b. Relative 
resolution of the nonunifonn DFT in Figure 3c is calculated as l/(2/ / ' / r s )=0.5 and it is half the nonnal 
resolution because of analysis is performed in two Nyquist zones. Nevertheless, squares under blue 
and black plots in Figure 3c are equal to one's depicted in Figure 2d. The maximum increase in the 
frequency resolution 2000/64*31 times is achieved on a complex exponent at frequency 0.35 Hz 
by the EDFT. 

(a) 



(b) 



(c) 



Figure 3. The estimates obtained in the extended frequency range: 

(a) Power spectrum of unifonn sequence - True (red), DFT (blue) and EDFT (black), 

(b) Power spectrum of nonuniform sequence - True (red), DFT (blue) and EDFT (black), 
(c) Relative frequency resolution of nonunifonn sequence - DFT (blue) and EDFT (black). 
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The EDFT also increases resolution in half to process transient and random signal components with 
the normal frequency resolution equal to 1, as it is indicated by the red dotted lines in Figure 3c. 
Hence the conclusion that EDFT can handle nonunifonnly sampled signals in multiple Nyquist 
zones only if the overall spectrum of band-limited signal components does not exceed one zone. 
Since the spectrum of unifonn sequence (red color line in Figure 1) does not cover the entire 
Nyquist zone the EDFT should be able to handle it with mean sampling period T s greater than T but 
less than 2 T. The increase of T s could be achieved by skipping samples from the unifonn sequence 
randomly. The resulting sequence is considered as nonunifonnly sampled because the distance 
between adjacent readings become unequal. The power spectra in Figure 4 show an example of the 
impact of sample skipping on the performance of DFT and EDFT. Input sequences are modeled by 
removing 16 and 24 samples randomly from the uniform 64-point data and leads to increase 
72=64/487"= 1,33s and 72=64/407= 1,6s, respectively. The simulation shows that DFT fails to 
process sequences with missing samples, while EDFT is still applicable (Fig.4.b) if one Nyquist 
zone limit on the total signal component spectrum is satisfied, otherwise the estimate becomes 
worse (Fig.4.c). Note that the result depends not only on the number of skipped samples, but also 
on their distribution within the sequence. The most sensitive to missing samples are transient signals 
which require dense sampling within their location, whereas deterministic signals appear more 
resistant, especially if the frequencies of discrete components lie on the EDFT grid. It is expected 
that the considerably greater increase of mean sampling period 72 can be achieved for pure 
deterministic signals [7]. 

(a) 





Figure 4. The power spectrum - True (red), DFT (blue) and EDFT (black), 
of the 64-point sequence without losses (a) and with randomly skipped (b) 16 and (c) 24 samples. 
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Let's validate above expectations for signal consisting only of sinusoids in a white noise and 
generate a real-value sequence of the length K=6 4 samples as a sum of four sine waves with 
predefined amplitudes 0.5, 1, 2 and 3, having arbitrary initial phases and randomly selected 
frequencies on the EDFT grid (0.001 Hz). Moreover, the signal sampling is also perfonned on the 
grid with T= 1 second lag by choosing 64 time points within NT= 1000 second interval on a random 
basis and results in about a 16-fold increase of the average sampling period, T S =NT/K= 15,675s. 
Finally, a white Gaussian noise with SNR=20 dB is added and nonunifonn real-value sequence 
illustrated in Figure 5 a. 

(a) 



0 100 200 300 400 500 600 700 800 900 1000 


Sample number 
(b) 



(c) 


-1-1 

-1-1- 

lIu. AiL a it .• j .,. 

- 1 i- 1 -r 

i ilM.. il J .. .iluil jl Ji.JrA .jL.JL.liAa A .i 

UjL'M 




0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

Frequency 

(d) 



100 200 300 400 500 600 700 800 900 1000 

Sample number 


Figure 5. Real-value sequence (a) - 64 nonunifonn samples, SNR 20 dB, 

Amplitude spectrum (b) - True values (red cycles) and estimates by DFT (blue), EDFT (black), 
Relative frequency resolution (c) by DFT (blue line) and EDFT (black), 

Original (blue cycles) and interpolated sequence (d) by Inverse EDFT - 1000 unifonn samples. 
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The true frequencies and amplitudes (red cycles) as well as the amplitude estimates of DFT (blue 
line) and EDFT (back line) are depicted in the Figure 5b and showed that DFT cannot recognize 
weaker power sinusoids while the Extended DFT picks up all of them and estimates their 
amplitudes and phases precisely. The performance difference is explained in the Figure 5c by 
comparison of the resolution of both DFTs with respect to the normal frequency resolution (equal 
to one). The relative resolution of the DFT (blue line) is calculated as 1 l{2f u T s )=KI N=D .D6A and it 
is considerably less than it is required for successful signal processing. This causes aliasing and 
leakage effects, because the spectrum of the sequence spreading in almost 16 Nyquist zones and N- 

K samples at the input of the DFT could be considered as zeroed by the rectangular windows. The 

1 1 

relative resolution of the EDFT (black line) is calculated as - F./S = - F./S and it increases 

2 fu T sK N 

NIK times reaching the value close to one at frequencies of sinusoids. Thus, signal processing with 
just a normal frequency resolution allows EDFT not only estimate the parameters of the signal 
components correctly, but also IDFT applied to its output returns a sequence of length N consisting 
of the original K and N-K interpolated samples (see Figure 5d). It should be noted that only a 
deterministic part of the signal is interpolated by EDFT, whereas a white Gaussian noise stays 
localized in time around the sampling points (blue cycles). 

The next sequence used in the computer simulations is well-known Marple&Kay data set taken 
from [3]. It is 64-point real sample sequence of a process consisting of two unit power sine waves 
with frequencies of 0.2 and 0.21 Hz, a third one with a power of 0.1 (20 dB down) at 0.1 Hz and a 

(a) 



(b) 



(c) 



Figure 6. The power spectrum obtained for Marple&Kay data set by 
(a) DFT, (b) EDFT, (c) HRDFT. 
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colored noise in the frequency range [0.2...0.5] Hz (see red color lines in Figure 6). The signal 
upper frequency is/,=0.5 Hz and the length of the DFT is selected N= 1000. Only 500 positive 
frequencies are displayed, because Marple&Kay sequence is a real-valued and negative 
frequencies, if they are depicted, gives a symmetrical pattern to zero frequency. The Figure 5 shows 
the power spectra of the DFT, EDFT and HRDFT approaches in a single picture, while separately, 
these plots have been presented in [5] and [9], The performance of other well-known spectral 
analysis methods for Marple&Kay data set could be found in [3], including Minimum Variance 
approach, named in the Section 5.1 as a traditional Capon filter (37). 

The simulation results in the Figure 6a,b demonstrate that the classical DFT and EDFT can evaluate 
not only the spectrum of sinusoids, but also the shape of continuous spectrum of other signal 
components, whereas HRDFT on Figure 6c is suitable mostly for the estimation of a line spectrum. 
The plot in Figure 6a shows that due to limited frequency resolution the classical DFT cannot 
resolve sine waves at the frequencies 0.2 and 0.21. Although the first EDFT iteration coincides with 
DFT, in the further iterations EDFT is able to increase the frequency resolution around the powerful 
signal components and all three sine waves are clearly distinguished after the 15 th iteration in the 
Figure 6b. 

All the three DFTs have one common feature - the ability to get back 64 samples of Marple&Kay 
data set by applying IDFT to the output of each of these methods. Since the length of the DFT is 
chosen equal to 1000, the inverse transfonn (29) returns 1000-64 additional samples, which are 

(a) 



Sample number 

(c) 



0 100 200 300 400 500 600 700 800 900 1000 

Sample number 


Figure 7. Marple&Kay sequence (blue) and extrapolated data (black) 
by inverse (a) DFT, (b) EDFT, (c) HRDFT. 
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plotted in Figure 7 (black). The samples 65, 66, 67, ... are considered as a forward extrapolation, 
but samples 1000, 999, 998,... as a backward extrapolation of known 64-sample sequence (blue). 
Of course, the Marple&Kay sequence outside of giving data set is unknown, and plots on Figure 6 
are just three possible versions of its extrapolation. The classical DFT (Fig.7a) suggests that 
Marple&Kay sequence outside of given 64 samples will be zeros, HRDFT (Fig.7c) shows that the 
extrapolated data even will increase in power, while EDFT (Fig.7b) expects that the sequence 
beyond will have approximately the same power, which only gradually decreases in time. 

(a) 
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(b) 
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(c) 
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Sample number 

Figure 8. White Gaussian noise (blue) and extrapolated data (black) 
by inverse (a) DFT, (b) EDFT, (c) HRDFT. 

At the end of the computer modeling, we check extrapolated sequences obtained at the output of 
IDFT if Marple&Kay input data are replaced by the same size white Gaussian noise (Figure 8). 
According to the theory the PSD of white Gaussian noise should be constant (flat) across the entire 
frequency range and the readings in a such sequence are uncorrelated random variables, therefore 
they cannot be extrapolated. In practice, because of finite length sequences and pseudo-random 
generators used in the simulations, the above expectations are satisfied only approximately. The 
classical DFT, as the case Marple&Kay data illustrated in Figure 7a, yields zeros outside of given 
64-point sequence also in Figure 8a, that this time is perfectly consistent with the theory. 
Extrapolate by the EDFT (Fig.8b) vanish quickly, and this still agrees with the theory if practical 
considerations are taken into account. The HRDFT (Fig.8c) in contrary to DFT and EDFT extends 
the white Gaussian noise up to a length of 1000 samples showing a strong correlation in the input 
sequence and this is very unlikely to be true. 
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Any approach that claims that it is a high frequency resolution method in accordance with the 
Uncertainty Principle must make certain assumptions about the data outside of the observation 
period even if by itself it is not able to recover the signal. The advantage of the proposed method 
over similar ones is that EDFT based on a solution that satisfies the minimum least squares criteria 
(6), making it an accurate, reliable and stable. 

Run MATLAB program EDFTFIG.m available on mathworks.com to recreate the computer 
simulations presented in this section. 

7 EDFT in MATLAB code 

The EDFT package consisting of programs written in a simple MATLAB code and created to 
demonstrate the Extended DFT capabilities described in the previous sections. Each function 
contains commented (%) help text section where its syntax, algorithm, usage and features are 
described. 

The programs NEDFT.m and the inverse transform INEDFT.m can be applied for unifonn or 
nonuniform input/output data and frequency sets. 

function [F,S,Stopit]=nedft(X,tk,fn,I,W) 

% NEDFT Nonuniform Extended Discrete Fourier Transform. 

% 

% SYNTAX 

% a. Mandatory inputs/outputs 
% F=nedft(X,tk,fn); 

% Function NEDFT returns discrete Fourier transform F(fn) of input sequence X sampled at arbitrary selected time moments tk, where 

% frequencies fn also may be selected arbitrary. If size of fn is less than size of X, inputs X and tk will be truncated. 

% b. Mandatory and optional inputs/outputs 
% [F,S,Stopit]=nedft(X,tk,fn,I,W); 

% I Optional input parameter I can be used for limiting maximum number of iterations. If I is not specified in input arguments, default 
% value for I is set by parameter 'Miteration', that is, nedft(X,tk,fn)=nedft(X,tk,fn,Miteration). To complete iterations faster, the value 

% for 'Miteration' should be decreased. 

% W Input weight vector W, if specified, override the default values W=ones(size(fn)). W has at least one nonzero element. 

% S The second output argument S represents the Amplitude spectrum. Peak values of abs(S) can be used for estimate amplitudes of 

% sinusoids in the input sequence X. 

% Stopit is an informative output parameter. The first row of 'Stopit' showing the number of performed iteration, the second row indicate 
% breaking of iteration reason and may have the following values: 

% 0 - Maximum number of iteration performed. 

% 1 - Sum of outputs division sum(F./S) is not equal to length(fn)*length(X) within Relative deviation 'Rdeviat'. The calculations were 

% interrupted because of results could be inaccurate. 

% 2 - Relative threshold 'Rthresh' reached. To complete iteration process faster, the value for 'Rthresh' should be increased. 

% ALGORITHM 
% Input: 

% X - input sequence; 

% E - complex exponents matrix (Fourier transform basis) - E=exp(-i*2*pi*tk.'*fn); 

% I - (optional) number of maximum iteration. 

% W - (optional) weight vector W. The default values W = ones(size(fn)). 

% Output F and S for each NEDFT iteration are calculated by following formulas: 

% 1. R=E*diag(W/length(fn))*E'; 

% 2. F=W. * (X * in v(R) * E); S=(X*inv(R)*E)Vdiag(E'*inv(R)*E).'; 

% 3. W=S.*conj(S); - the weight vector W for the next iteration. 

% If length(X) is equal to length(fn), the NEDFT output do not depend on selected weight vector W and is calculated in non-iterative way. 

% TIPS for selection of mandatory NEDFT inputs X(tk) and fn: 

% 1. Input sequence X(tk) for NEDFT can be sampled uniformly or nonuniformly. Uniform sampling can be considered as a special 

% case of nonuniform sampling, where tk=[0,l,...,K-l]*T and T is sampling period. Nonuniform sampling can be realized in many 

% different ways, like as: 

% - uniform sampling with randomly missed samples (known as sparse data); 

% - uniform sampling with missed data segments (known as gapped data); 

% - uniform sampling with jitter: tk=([0,l,...,K-l] + jitter*rand(l,K))*Ts, where value for jitter is selected in range [0... 1 [ and Ts is 

% the mean sampling period; 

% - additive nonuniform sampling: tk=tk-l + (l+jitter*(rand-0.5))*Ts, k=l,...,K-l, t0=0; 

% - signal dependent sampling, e.g, level-crossing sampling, etc.... 

% 2. Frequencies for fn can be selected arbitrary. This mean, that user can choose not only the length of NEDFT (number of frequencies 

% in fn), but also the way how to distribute frequencies along the frequency axis. On other hand, to get adequate sequence X 

% representation, frequencies fn should be selected to cover overall range, where the input sequence X spectrum is supposed to be 

% found, otherwise, in result of NEDFT, all components having spectra outside fn will be incorporated. Note that fn should contain 

% negative frequencies too, and for a real value X(tk) analysis each positive frequency in fn should have corresponding negative one, 

% for example fn=[-ceil((N-l)/2):floor((N-l)/2)]/(N*T), where N=length(fn). 

% 3. Frequencies for vector fn can be added in any order. Therefore it is possible to combine different frequency sets in one or just add 

% individual frequencies of interest to fn for example, fn=[fnl fn2 fl f2], where fnl and fn2 are different frequency sets, fl ,f2 - specific 
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F(fn)=[F(fnl) F(fn2) F(fl) F(£2)], S(fn)=[S(fnl) S(fn2) S(fl) S(f2)]. 


% frequencies. NEDFT outputs will be calculated accordingly ■ 

% FEATURES 

% 1. NEDFT output F(fn) is the discrete Fourier transform of sequence X(tk). The Power Spectral Density function of nonuniform 

% sequence X(tk) can be estimated by the following fonnula: abs(F). A 2/(N*Ts). 

% 2. In general, the function Y=inedft(F,fn,tn) (see attached program) is used to calculate the reconstructed sequence Y(tn). 

% If frequencies fn are selected on the same grid as used by FFT algorithm, then ifft(F) can be applied to get uniformly re-sampled and 

% extrapolated to length(fn) version of input sequence X(tk). 

% 3. NEDFT output S(fn) estimate amplitudes and phases of sinusoidal components in sequence X(tk). 

% 4. NEDFT can increase frequency resolution length(fn)/length(X) times. Division of outputs l/(Ts*(F./S)) demonstrate the frequency 

% resolution of NEDFT. The following is true for any NEDFT iteration: 0<F./S<=length(fn) and sum(F./S)=length(fh)*length(X). 

% 5. If input arguments are matrixes, the NEDFT operation is applied to each column. 

% See also INEDFT, IFFT, EDFT. 

% Author: Vilnis Liepins (vilnislp@gmail.com) 

% Reference: V.Liepins, "A spectral estimation method of nonuniformly sampled band-limited signals. Automatic Control and Computer 
% Sciences,Vol.28, No.2, pp.66-77, 1994. 

% 

% ~ --- Set default parameters for NEDFT = - = - == - ===== 


% Limit for maximum number of iteration (Stopit 0). 
% Value for relative deviation (Stopit 1). 

% Value for relative threshold (Stopit 2). 

= Check NEDFT input arguments ==== =-=-=-= == 


Miteration=30; 

Rdeviat=0.0005; 

Rthresh=0.0001; 

% .. 

if nargin<3,error('Not enough input arguments. See help nedft.'),end 
if sum(any(isinf(X)))|sum(any(isnan(X))), error('Input argument X contain Inf or NaN. See help nedft.'),end 
if size(X,l)==l,trf=0;else, X=X.'; tk=tk.'; fh=fn.'; trf=l;end % Check input argument X. 

[L K]=size(X); % K - the length of X. 

if size(tk,l)~=L | size(tk,2)~=K, error('Size of input arguments X and tk must be equal. See help nedft.'),end 
if size(fn,l)~=L, error('Incorrect size of input argument fn. See Help nedft.'),end 
N=size(fn,2); % N - the length of DFT. 

if N<K,X=X(:,1:N); tk=tk(:,l:N); K=N;end % Truncate sequence X if N<K. 

if nargin<4,I=Miteration;else, % Set value for maximum number of iterations, 

if isempty(I),I=Miteration;end,I=floor(I(l));end % Check input argument I. 

if nargin>4,if trf==l ,W=W.';end % Check input argument W. 

if (size(W,2)~=N)|(size(W,l)~=L),error('Incorrect size of input argument W. See help nedft.'),end 
W=W.*conj(W); % W have elements >=0. 

if sum(any(isinf(W)))|sum(any(isnan(W))), error('Input argument W contain Inf or NaN. See help nedft.'),end 
if any(find(sum(W>0)<l)), error('Too many zeros in input argument W. See help edft.'fend 


else,W=ones(L,N);end 

%====== 


% Default values for W. 


W=ones(L,K); end 


if K=N, 1=1; 

% .. = 

F=zeros(L,N); S=zeros(L,N);it=l; 
Stopit=[I*ones(l ,L);zeros(l ,L)]; 

%====== 


== Check for a special case ===================== 

% If K=N, perform just one NEDFT iteration. 
Set default values for NEDFT output arguments ========= 

% Fill zeros in output matrixes F and S. 

% Stopit 0: Set values for default Stopit. 

=== Calculate NEDFT for each X column 1 = ===================== 


% Calculate the complex exponents matrix E. 
% Start iterations... 


for 1=1 :L, 

E=exp(-i*2*pi*tk(l,:).'*fn(l,:)); 
for it=l :I, 

% Calculate the correlation matrix R by using vectorized form or 
% R=E*diag(W(l,:)/N)*E'; 

for n=l :K, % a loop structure 

for k=n:K, 

R(k,n)=sum(W(l,:). *conj (E(n,:)). *E(k, :))/N; 
if n=k,R(n,k)=conj(R(k,n));else,R(n,n)=real(R(n,n));end 
end 
end 

% Calculate RE=R\E or pinv(R)*E if R is badly conditioned. 

if rcond(R)>eps, RE=R\E;else,RE=pinv(R)*E;end 
% Calculate ERE=diag(E'*inv(R)*E).'=sum(conj(E).*RE). 

ERE=sum(conj(E).*RE); 

% Calculate outputs for iteration (it): N-point NEDFT (F) and Amplitude Spectrum (S). 
F(1,:)=X(1,:)*RE; 

S(1,:)=F(1,:)./ERE; 

F(1,:)=F(1,:).*W(1,:); 

% Stopit 1: Break iterations if sum(F./S) is not equal to N*K. 

if abs(ERE*W(l,:).'/N/K-l)>Rdeviat,Stopit(:,l)=[it;l]: 
if sum(W(l,:)=0)>0, S(l,W(l,:)=0)=zeros(l,sum(W(l,:)==0));end,break,end 
% Calculate weight (W) for the next iteration. 

W(l,:)=S(l,:).*conj(S(l,:)); 

% Stopit 2: Break iterations if relative threshold reached. 

SW(it)=sum(W(l,:)); 
if it>l ,thit=abs(SW(it-1 )-SW(it))/SW(l); 
if thit<=Rthresh,Stopit(:,l)=[it;2];break,end,end 
end % ... end iterations. 


end 

% .= - 

if trf=l,F=F.';S=S.';end 


= Adjust size of NEDFT output 
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function Y=inedft(F,fn,tn) 

% INEDFT Inverse Nonunifonn Extended Discrete Fourier Transform. 

% 

% Y=inedft(F,fn,tn) is the inverse discrete Fourier transfonn of the vector F estimated by NEDFT function at arbitrary frequency 

% set fh: F(fn) » Y(tn), where time moments tn for reconstructed sequence Y can be uniformly or nonuniformly spaced in time. 

% In the special case of uniform vectors fh and tn, the INEDFT function can be replaced by well-known MATLAB function IFFT. 

% If input arguments are matrixes, the INEDFT operation is applied to each column. 

% 

% See also IFFT, EDFT, NEDFT. 

% 

% Vilnis Liepins (vilnislp@gmail.com) 

% 

%== - === Check INEDFT input arguments ========================== 

if nargin<3,error('Not enough input arguments. See help inedft.'),end 
if size(F,l)=l, 
trf=l; F=F.'; tn=tn.'; 
else 

trf=0; fh=fh.'; 
end 

[N L]=size(F); 

if size(fn,2)=N, errorfSizes of input arguments F and fh must be equal. See help inedft.'),end 
if size(tn,2)~=L, error('Incorrect size of input argument tn. See help inedft.'),end 

%=== ■■—.-■■■■ ' ■.■ .. Calculate INEDFT for each X column 1 ======================= 

for 1=1 :L 

E=exp(i*2*pi*tn(:,l)*fn(l,:)); 

Y(:,1)=E*F(:,1)/N; 

end 

% - = - ===== - Adjust size of INEDFT output ================================ - 

if trf=l ,Y=Y.';end 


From the viewpoint of calculation complexity is reasonable to use the same frequency grid as 
Fast Fourier Transfonn (FFT.m in MATLAB library). This allows to apply the FFT algorithm 
in EDFT calculations, which considerably reduce computation time, because each FFT requiring 
number of operations proportional to Mog(A) rather than N 2 [2], however the efficiency of FFT 
algorithm could be in between of these two values since it depends also on the value of N. 
EDFT.m program is designed as a faster realization of Extended DFT, where the algorithm 
described in [5,7] is implemented. The code is applicable for uniformly sampled signals or data 
with gaps in it. The inverse transfonn to EDFT.m is MATLAB library program IFFT.m. 

function [F,S,Stopit]=edft(X,N,I,W) 

% EDFT Extended Discrete Fourier Transform. 

% 

% Function EDFT produce discrete N-point Fourier transform F and amplitude spectrum S of the datavector X. Data X may contain NaN 
% (Not-a-Number). Inverse transform is Matlab built-in function ifft. 

% 

% SYNTAX 

% [F,S,Stopit]=edfi(X,N) for N>length(X) calculate F and S iteratively (see ALGORITHM below). If data X do not contain NaN and 
% N<=length(X)or N is not specified, EDFT return the same results as fast Fourier transform: F=fft(X,N) or F=ffi(X) and S=F/N. 

% [F,S,Stopit]=edfi(X,N,I) performs edft(X,N) with limit I for maximum number of iterations. Default value for I is set by parameter 
% 'Miteration', that is, edft(X,N)=edft(X,N,Miteration). To complete iteration process faster, the value for 'Miteration' should 

% be decreased. 

% [F,S,Stopit]=edfi(X,N,I,W) execute edft(X,N,I) with initial conditions defined by weight vector W. Default values for W are ones(size(F)). 
% W must have at least length(X) nonzero elements. 

% Stopit is an informative (optional) output parameter. The first row of Stopit showing the number of performed iteration, the second row 
% indicate breaking of iteration reason and may have the following values: 

% 0 - Maximum number of iteration performed. If length(X)<=N, only one EDFT iteration is perfonned (1=1). 

% 1 - Sum of outputs division sum(F./S) is not equal to length(X)*N within Relative deviation 'Rdeviat'. The calculations were 

% interrupted because of results could be inaccurate. If this occur in the first iteration, then outputs F and S are zeros. 

% 2 - Relative threshold 'Rthresh' reached. To complete iteration process faster, the value for 'Rthresh' should be increased. 

% ALGORITHM 
% Input: 

% X - input data. 

% N - length of discrete Fourier transfonn. 

% I - (optional) number of maximum iteration. If not specified, I=Miteration. 

% W - (optional) weight vector W. If not specified, W = ones( 1 ,N) is used for the first iteration. 

% E - Fourier transform basis matrix: E=exp(-i*2*pi*(0:length(X)-l)'*(0:N-l)/N); 

% If part of unknown data in X are replaced by NaN then the time vector (Odength(X)-l) is changed to exclude time moments where 

% NaNs inserted. 

% Output F and S for each EDFT iteration are calculated by following formulas: 

% 1. R=E*diag(W/N)*E'; 
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% EDFT using function ifft to calculate R faster. 

% 2. F=W.*(X*inv(R)*E); S=(X*inv(R)*E)./diag(E'*inv(R)*E).'; 

% Levinson-Durbin recursion used to inverse toeplitz R. Function fit applied to speed up matrix multiplications. 

% 3. W=S.*conj(S); W used as input to the next EDFT iteration. 

% A special case: If length(X) is equal to N, the EDFT output do not depend on selected weight vector W and is calculated in non-iterative way. 
% FEATURES 

% 1. EDFT output F is the N-point Fourier transform of data X. The Power Spectral Density (PSD) function can be calculated by the 

% following formula: abs(F). A 2/(N*T), T - sampling period. 

% 2. EDFT can extrapolate input data X to length N. That is, if apply EDFT for N>length(X), get the results: 

% F=edft(X,N)=edft(Y)=fft(Y); Y=ifft(F);, where Y is input X plus non-zero forward and backward extrapolation of X to length N. 

% 3. EDFT output S estimate amplitudes and phases of sinusoidal components in input data X. 

% 4. EDFT can increase frequency resolution N/length(X) times. Division of outputs 1/(T*F./S) demonstrate the frequency resolution 

% of EDFT. The following is true or any EDFT iteration: 0<F./S<=N, sum(F./S)=N*length(X). 

% 5. EDFT input data X may contain NaN. NaNs indicate unavailable data or missing samples or data segments in X. Outputs F and 

% S are calculated by applying slower algorithm then in case of X without NaN. 

% 6. If X is a matrix, the EDFT operation is applied to each column. 

% See also FFT, IFFT, FFTSHIFT. 

% 

% Vilnis Liepins (vilnislp@gmail.com) 

% Reference: V. Liepins, "An algorithm for evaluating a discrete Fourier transform for incomplete data", Automatic control and computer 
% sciences, Vol.30, No.3, pp.27-40, 1996. 

% NOTE: The first version of file (gdft.m) was submitted on 10/7/1997 as Matlab 4.1 code. 


= Set default parameters for EDFT ================== 

% Limit for maximum number of iteration (Stopit 0). 
% Value for relative deviation (Stopit 1). 

% Value for relative threshold (Stopit 2). 

= Check EDFT input arguments ============================ 


%====”= 

Miteration=30; 

Rdeviat=0.0005; 

Rthresh=0.0001; 

% .. = 

if nargin==0, error('Not enough input arguments. See help edft.'), end 
if sum(any(isinf(X))), error('Input argument X contain Inf. See help edft.'),end 
if size(X,l)==l,X=X.';trf=l;else, trf=0; end % Check input argument X 

[K L]=size(X); % K - length of input argument X 

if nargin>l,if isempty(N),N=K;end,N=floor(N(l));% Check input argumentN. 

if N<K, X=X( 1 :N,: );K=N;end,else,N=K;end % Truncate X if has more than N points 

Xnan=~isnan(X); % Check X on NaNs: 

if N=l, KK=Xnan; % Xnan - indicate data as T , NaN as '0' 

else, KK=sum(Xnan);end % KK - length of input X without NaN 

if nargin<3, I=Miteration; % Check input argument I. 

else,if isempty(I),I=Miteration;end,I=floor(I(l));end % Set default value for I. 

if nargin<4,W=ones(N,L);else,if trf=l ,W=W.';end % Set default values for W. 

if (size(W,l)~=N)|(size(W,2)~=L), error('Incorrect size of input argument W. See help edft. 1 ),end 

W=W.*conj(W); % W have elements >=0 

if any(find(sum(W>0)<KK)), error('Too many zeros in input argument W. See help edft.'),end,end 


%======= 

F=zeros(N,L);S=zeros(N,L); 
Stopit=[I*ones( 1 ,L);zeros( 1 ,L)]; 
% .. = 

for 1=1 :L, 

%============= 


= Set default values for EDFT output arguments ========== 

% Fill with zeros output matrixes F,S. 
% Set default value for Stopit. 

= Calculate EDFT for each X column 1 ================ 


= Check for a special cases = 


if KK(1)=N|KX(1)==0, 

% If length(X)=N or X(:,l) has all NaNs then 

F(:,l)=fft(X(:,l),N); 

% EDFT output (F,S) equals to FFT. 

S(:,1)=F(:,1)/N; 


Stopit(:,l)=[l; 0]; 


elseif K=1&N~=1, 

% Special case, the length(X)=l, 

F(:,l)=ffi(X(:,l),N).'; 

% EDFT output (F,S) equals to FFT. 

S(:,1)=F(:,1)/N; 


Stopit(:,l)=[l; 0]; 


elseif isempty(find(X(:,l)))&KK(l)>0, 

% If input X(:,l) has all zeros or zeros&NaN 

Stopit(:,!)=[!; 0]; 

% then EDFT output (F,S) is zeros. 


%==== 

elseif KK(1)==K, 

%=== 


= Basic EDFT algorithm started = 


% hiput X(:,l) does not contain NaN 


========= Apply FASTER algorithm = 

for it=l :l, % Start iterations... 

% Calculate correlation vector (r). 
r=ifft(W(:,l)); 

% Perform inverse of correlation matrix R by applying Matlab \ operator or 
% a=[l; toeplitz(conj(r(l:K-l)))\(-r(2:K))];V=a.'*conj(r(l:K)); 

a=-r(2)/r(l); V=r(l)-r(2)*conj(r(2))/r(l); % Levinson-Durbin recursion, 

for n=l:K-2, 

alfa=[l a.']*r(n+2:-l:2);rho=-alfa/V;V=V+rho*conj(alfa);a=[a+rho*conj(flipud(a));rho]; 
end 


a=[l;a]; 

% Calculate ERE=diag(E'*inv(R)*E) and XR=X*inv(R). 
XR=zeros(K, 1 );RE=zeros(K, 1 );rc=a; 
for k=l:K/2, 
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kO=K-k+l ;kl=2:K-2*k+l; k2=k+l:K-k; k3=k:K-k+l; 

RE(l)=RE(l)+2*rc(k);RE(kO-k+l)=RE(kO-k+l)+2*rc(kO);RE(kl)=RE(kl)+4*rc(k2); 

XR(k)=XR(k)+rc(k3) , *X(k3,l);XR(kO)=XR(kO)+(flipud(rc(k3))).'*X(k3,l); 

XR(k2)=XR(k2)+rc(k2)*X(k,l)+flipud(conj(rc(k2)))*X(kO,l); 

rc(k2)=rc(k2-l)+conj(a(k+l ))*a(k2)-a(k0)*flipud(conj(a(k2+l))); 

end 

if round(K/2)>K/2, 

RE( 1 )=RE( 1 )+rc(k+1 );XR(k+l )=XR(k+1 )+X(k+1 ,l)*rc(k+1); 
end 

ERE=real(fft(RE,N));W(:,l)=W(:,l)/real(V); 

% Stopit 1: Break iterations if sum(F./S) is not equal to N*K or NaN. 

stit=abs(ERE.'*W(:,l)/N/K-l);if (stit>Rdeviat)|isnan(stit), Stopit(:,l)=[it-l; 1]; break,end 
% Calculate outputs for iteration (it): N-point EDFT (F) and Amplitude Spectrum (S). 
F(:,l)=fft(XR,N); 

S(:,1)=F(:,1)./ERE; 

F(:,1)=F(:,1).*W(:,1); 

% Calculate weight (W) for the next iteration. 

W(:,l)=S(:,l).*conj(S(:,l)); 

% Stopit 2: Break iterations if relative threshold reached. 

SW(it)=sum(W(:,l)); 

if it> 1, thit=abs(SW(it-l)-SW(it))/SW(l);if thit<=Rthresh, Stopit(:,l)=[it; 2]; break,end 
end 

end % ... end iterations. 


%== -. — .. - : E n d 0 f FASTER algorithm ====================== 

else % Input X(:,l) contains NaN 

%== -~-== = = = Apply SLOWER algorithm ======================== 

INVR=zeros(K);ER=zeros(K, 1); 

X(find(~Xnan(:,l)),l)=zeros(K-KK(l),l); % Replace NaN by 0 in X 

t=fmd(Xnan(:,l)); % Sample numbers vector (t) 

for it=l :I, % Start iterations... 

% Calculate correlation matrix R by applying ifft and inverse R by Matlab backslash operator. 
RT=iffi(W(:,l)); 

R=toeplitz(RT( 1 :K)); 

INVR(t,t)=R(t,t)\eye(KK(l)); 

% Calculate ERE=diag(E'*inv(R)*E).' by applying fft. 

ER( 1 )=trace(INVR); 
for k=l:K-l 

ER(k+l, 1 )=sum(diag(INVR,k)+conj(diag(INVR,-k))); 
end 

ERE=real(ffl(ER,N)); 

% Stopit 1: Break iterations if sum(F./S) is not equal to N*KK. or NaN. 

stit=abs(ERE.'*W(:,l)/N/KK(l)-l);if (stit>Rdeviat)|isnan(stit), Stopit(:,l)=[it-l; l];break,end 
% Calculate outputs for iteration (it): N-point EDFT (F) and Amplitude Spectrum (S). 
F(:,l)=fft(conj(INVR)*X(:,l),N); 

S(:,1)=F(:,1)./ERE; 

F(:,1)=F(:,1).*W(:,1); 

% Calculate weight (W) for the next iteration. 

W(:,l)=S(:,l).*conj(S(:,l)); 

% Stopit 2: Break iterations if relative threshold reached. 

SW(it)=sum(W(:,l)); 

if it>l, thit=abs(SW(it-l)-SW(it))/SW(l);if thit<=Rthresh,Stopit(:,l)=[it; 2];break,end 
end 

end % ... end iterations. 

%== .~ . == = = = End of SLOWER algorithm ======================= 

end 

end 


% - = - - == - Adjust size of EDFT output 

if trf=l,F=F.';S=S.';end 


The next program demonstrates the applicability of the Extended DFT in 2-dimensional signal 
processing. The EDFT2.m program is based on the MATLAB library program FFT2.m where 
FFT.m calls are simply replaced by EDFT.m. The inverse transform to EDFT2.m is the 
MATLAB library program IFFT2.m. 

function f = EDFT2(x, mrows, ncols) 

% EDFT2 Two-dimensional Extended Discrete Fourier Transfonn. 

% 

% EDFT2(X) returns the two-dimensional Fourier transfonn of matrix X. Before running EDFT2 unknown data (if any) inside of X 

% should be replaced by NaN (Not-a-Number). If X is a vector, the result will have the same orientation. 

% EDFT2(X,MROWS,NCOLS) performing size MROWS-by-NCOLS Fourier transform without padding of matrix X with zeros. 

% Inverse transfonn is Matlab built-in function ifft2. 

% See also EDFT, IFFT2. 

% 
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% EDFT2 is created on basis of FFT2 (J.N. Little 12/18/1985) by Vilnis Liepins. 
% No input. 

if nargin==0, error('Not enough input arguments. See help edft2.'), end 
[m, n] = size(x); 

% Basic algorithm. 

if (nargin == 1) & (m > 1) & (n > 1) 

% f = fft(fft(x).').'; 
f = edft(edft(x).').'; 
return; 
end 

% Padding for vector input. 

if nargin < 3, ncols = n; end 

if nargin < 2, mrows = m; end 

mpad = mrows; npad = ncols; 

if m = 1 & mpad > m, x(2, 1) = 0; m = 2; end 

if n == 1 & npad > n, x(l, 2) = 0; n = 2; end 

if m = 1, mpad = npad; npad = 1; end % For row vector. 

% Transform. 

%f= fft(x, mpad); 

%if m > 1 & n > 1, f = fft(f.', npad).'; end 
f = edft(x, mpad); 

ifm> 1 & n> l,f = edft(f.', npad).'; end 


The first version of EDFT (file GDFT.m) was submitted to file-exchange server on 10/7/1997 
as MATFAB 4.1 code. The renewed code version uploaded on 8/5/2006 and available online 

mathworks.com and researchgate.net . 

Please note that programs have not been tested on the latest MATFAB versions and therefore 
have opportunities to performance improvements (see for example [ 25 , 28 ]). 
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